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1 Introduction 

There has been a resurgence of interest in charmonium spectroscopy over the last decade 
driven by experimental observations at the i?-factories, CLEO-c, the Tevatron and BES of 
rather narrow states close to or above open-charm thresholds. Some had been predicted but 
hitherto unobserved whereas many were unexpected, for example the enigmatic X(3872), 
y(4260) and other "X,Y,Z"s. Prior to these discoveries, the relatively small number of 
observed mesons all fitted into the simple pattern predicted by quark models. The new 
states have spurred theoretical discussion as to their nature with suggestions including 
conventional quark-antiquark states appearing at unexpected masses, hybrid mesons in 
which the gluonic field is excited, molecular mesons and tetraquarks. To date, however, 
there has been no observation of a state in the charmonium region with manifestly exotic 
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J {J, P and C are respectively spin, parity and charge-conjugation quantum numbers); 
an observation of such an exotic would be a smoking gun for physics beyond quark potential 
models. Experimental interest continues with BESIII, experiments at the LHC, and the 
planned PANDA experiment at GSI/FAIR. An extensive review of the experimental and 
theoretical situation is given in Ref. [1]. 

From this perspective, the charmonium system is a particularly interesting probe of 
the strong regime of QCD. Being reasonably non-relativistic it is well described below the 
open-charm {DD) threshold by modelling it as a quark and anti-quark bound together in 
a confining potential. It has also been studied using effective field theory approaches and 
in the framework of QCD sum rules. No real consensus has emerged about the nature of 
many of the new states. 

A vital part of understanding and testing QCD is to make precise predictions of the 
hadron spectrum it implies, as opposed to the spectrum in a QCD-inspired approach, and 
to test these predictions against high-quality experimental data. The charmonium system 
provides an important arena for this and, in particular, discerning the puzzling states within 
QCD is a crucial aspect of understanding the strong interaction. Lattice QCD provides 
a method for performing such ab-initio calculations within QCD: the theory is discretised 
on a finite lattice and then energies and other properties are extracted from numerical 
computations of Euclidean-time correlation functions. 

The precision calculation of the masses of the lower-lying states has long been an 
important benchmark of lattice calculations; examples of recent studies of the lightest 
charmonium mesons are given Refs. [2-5] and reviewed in Ref. [6]. More recently, progress 
has been made in using lattice QCD to study higher- lying states. A pioneering quenched 
calculation of the excited charmonium spectrum was presented in Ref. [7] and there have 
since been some dynamical studies [8, 9]; we comment on these previous studies later when 
we discuss the results of this work. 

In lattice calculations, the theory is discretised on a finite hypercubic grid and the 
full rotational symmetry of the continuum is reduced to the symmetry group of a cube. 
States at rest are then not classified according to spin, J, but by the irreducible represen- 
tations, irreps, of this reduced symmetry group; the various components of high spin states 
are distributed across several irreps. The Hadron Spectrum Collaboration has overcome 
the challenges this reduced symmetry presents by using a large basis of carefully con- 
structed operators with a combination of techniques and dynamical anisotropic lattices. 
This methodology has been shown to be very effective at extracting extensive spectra of 
light isovector and isoscalar mesons [10-12] and baryons [13-15], including highly excited 
states, states with high spin and those with exotic quantum numbers, and in reliably 
determining the continuum J^*-^ of the extracted states. 

In this article we present a dynamical lattice calculation of the excited spectrum of 
charmonium mesons, the first application to charmonium of these techniques. The results 
presented here go beyond any previous dynamical calculation of the excited charmonium 
spectrum. We will show how we can reliably extract an unprecedented number of states, 
including the spectrum of exotics below around 4.5 GeV, and identify their continuum 
quantum numbers. In particular, identification of the J^^ will be vital in order to com- 
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pare with experiment and so to disentangle the nature of the many puzzhng resonances. 
Prehminary results from this work have been presented in Ref. [16]. 

2 The anisotropic lattice action 

This study aims to determine the spectrum of charmonium states, including excitations 
and any states with an intrinsic gluonic component. In Euclidean space-time, this spectrum 
can be computed by observing the behaviour of correlation functions of appropriately con- 
structed operators, which fall exponentially as the time separation is increased with rate 
in proportion to the energies of the eigenstates of the Hamiltonian. Since the charmo- 
nium spectrum lies above 3 GeV, these correlation functions decrease rapidly at large 
Euclidean time separations. Correlation functions of creation operators which excite high- 
lying states also suffer from substantial statistical variance in Monte Carlo determinations. 
The anisotropic lattice provides a very useful framework to ameliorate these technical prob- 
lems, as it enables the correlator to be determined with a finer temporal resolution at a 
manageable computational cost. 

The anisotropic lattice offers a further advantage when studying charmonium physics. 
The charm quark mass, rric, has proved problematic in simulations in the past. It is too light 
for non-relativistic approaches to be valid whilst it is also large enough that aruc <C 1 only 
for the finest isotropic ensembles currently available. A number of successful approaches 
have been developed, including the Fermilab and HISQ actions [2, 17]. For this study, 
we will use a discretisation in which the spatial lattice spacing, a^, and temporal lattice 
spacing, at, are related via ^ = ag/at ~ 3.5. The fine temporal discretisation will ensure 
that we carry out our calculation with atnic ^ 1 and the standard relativistic formulation 
of fermion actions can be used to study states where the typical spatial momentum of a 
charm quark inside a hadron at rest is small. This should be a reasonable approximation 
for the spectroscopy measurements we make in the current work. Nevertheless, in the 
simulations carried out here, a^mc is also less than unity. We explore the size of 0{as) 
errors in Section 5.3. 

2.1 The lattice action 

For this study, we make use of the substantial ensembles of anisotropic gauge-field con- 
figurations generated by the Hadron Spectrum Collaboration [18, 19]. The gauge ac- 
tion is tree-level Symanzik-improved while in the fermion sector we use the anisotropic 
Shekholeslami-Wohlert action with tree-level tadpole improvement and three-dimensional 
stout-link smearing [20] of gauge fields. The fields are drawn from the importance sam- 
pling distribution appropriate for studies with two dynamical light quarks and a dynamical 
strange quark (A'^^ = 2 + 1); more details are given in Refs. [18, 19]. 

The charm quark action is the same as that used for the light and strange quarks, 
except that the parameter describing the weighting of spatial and temporal derivatives 
in the lattice action is chosen to give a relativistic dispersion relation for low-momentum 
mesons containing charm quarks. This non-perturbative determination [18] is described 
in more detail in Section 2.2. The charm quark fields are not treated dynamically in the 



-3- 



Volume m^/MeV iVcfgs Nt 



tsrcs 



vocs 



16^ X 128 396 96 
24^ X 128 396 552 



128 



32 



162 



64 



Table 1. The gauge-field ensembles used in this work. A'cfgs and A'tsrcs are respectively the number 
of gauge-field configurations and time-sources per configuration used; N^ecs refers to the number of 
eigenvectors used in the distillation method [22] . 

configuration generation stage and so the quenching of these fields introduces a systematic 
uncertainty in our analysis. We expect the efi^ects on the spectrum of this non-unitary 
treatment of charm quarks to be fairly small. The bare mass parameter for the valence 
charm quarks is determined by ensuring the ratio of the r/c meson and the Q baryon masses 
takes its experimental value. 

To quote energies and lengths in physical units we follow Ref. [19] and consider the 
ratio of the mass of the baryon measured on these ensembles, atmn = 0.2951(22) [14], 
to the experimental mass, nin = 1672.45(29) MeV [21]. Setting the scale in this way, we 
find the parameters in the lattice action used in this study correspond to a spatial lattice 
spacing of ~ 0.12 fm and a temporal lattice spacing approximately 3.5 times smaller, 
a^^ ~ 5.7 GeV. The lattices employed here have m.,^ 400 MeV. Two volumes are used, 
{L/asY ^ (^/'^t) = ^ 128 and 24^ x 128, corresponding to spatial extents of ~ 1.9 fm 
and ~ 2.9 fm respectively. The ensembles used in this study are summarised in Table 1 
and full details of the lattice actions and parameters are given in Refs. [18, 19]. 

2.2 Determining the action parameters for the charm quark 

On the anisotropic lattice, the weights of any action terms that represent temporal and 
spatial derivatives must be determined in order that a physical observable takes its exper- 
imental value when calculated in the simulation. This is of particular importance for dy- 
namical simulations to ensure the correct continuum limit [23]. For these ensembles where 
the target ratio of scales is ^ = 3.5, this determination for the gauge fields along with the 
dynamical light and strange quarks is presented in Ref. [18]. This ratio has been measured 
in the light meson sector from the pion dispersion relation, giving = 3.444(6) [24]. 

As part of this study, we determined the dispersion relation for the r]c meson on the 
16^ ensemble. The dispersion relation for meson A is 



and since all quark fields making up the meson obey periodic boundary conditions for 
the finite spatial extent of the lattice, the momentum values are quantized, so asP = 
^{ux, riy, Hz) with m G {0, 1, . . . , L/us — 1}. Since we plan to follow up this work with a 
study of the scattering of mesons including charm quarks and hope to investigate the physics 
of states above the open-charm threshold, it is important to check that the dispersion 
relation is correct and that ^ is consistent with that in the light quark sector. 




(2.1) 
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Figure 1. Left panel: the dispersion relation of the rjc meson measured on 96 configurations for 
two choices of the bare parameter, 7/, in the charm quark action. As described in the text, 7; is 
the parameter in the light sea-quark action and 7c is determined by ensuring the correct dispersion 
relation for the rjc meson. Right panel: the dispersion relation for the heavy-light D meson measured 
on 39 configurations with = 7c in the charm quark action. Both on 16'^ volumes. 



In the left panel of Fig. 1, the dependence of the r/c meson energy with increasing 
quantized momenta, labelled by = -|- -|- , is plotted, along with the best straight- 
line fits. The triangular symbols show the dispersion relation when the bare parameter, 
■jf, in the charm quark action is the same as that used for the light sea quarks, 7/ = 7/ = 
3.4 [18]. Prom this data, the measured ratio of scales is = 3.18(2), a discrepancy from 
the target value of about 9%. 

To correct for this discrepancy, simulations were performed to determine the charm- 
quark action parameters that recover a consistent dispersion relation. Increasing the bare 
parameter to 7/ = 7c = 3.988, gives the second dispersion relation denoted by squares in 
the left panel of Fig. 1. The fit now yields a measured value of ^rjc = 3.50(2), consistent 
with the target value. Note that the statistical uncertainty on this determination is below 
the percent-level. We employ this 7/ = 7c in the charm quark action in all subsequent 
computations. 

As mentioned above, our future plans include a study of scattering of D mesons, so 
it is important to check whether the self-consistent action parameters presented above 
also result in a relativistic dispersion relation for the D meson. The light quark action 
parameters are left unchanged from those in the sea. This dispersion relation is presented 
in the right panel of Fig. 1 and the resulting slope yields a determination of the ratio 
of scales as = 3.38(3), just 3% different from the target value. The D meson mass 
extracted from this analysis is = 1893(1) MeV, but note that the simulated light 
quark mass is greater than its physical value, itItz ~ 400 MeV. 

We also note that the data in each of the figures agree with the relativistic dispersion 
relation up to n = (2,0,0). All five momenta are used in the linear fits which have 
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3 Operator construction and correlator analysis 

In lattice calculations, meson masses are extracted from two-point correlation functions, 

a,{t) = (0|O,(t)O](0)|0) , (3.1) 

where the interpolating operators Oj(0) create the state of interest at i = and Oi{t) 
annihilate the state at a later Euclidean time t. Inserting a complete set of eigenstates of 
the Hamiltonian we obtain the spectral decomposition, 

^'^■(*) = E^^"''"*' (3-2) 
where the vacuum-state matrix elements, 

Zf^{n\Ol\0) , (3.3) 

are referred to as overlaps, and the sum is over a discrete set of states because the calculation 
is performed in a finite volume. It is essential that the operators overlap well with the states 
under consideration. A well-established way to achieve this is to apply a smearing to the 
quark fields in the operators to emphasise the relevant low-energy modes. 

We use the distillation technique, described in Ref. [22], which implements a smearing 
and provides an efficient method to calculate correlation functions with large bases of 
operators. Furthermore, it allows us to project onto definite momentum at both the source 
and the sink with reasonable computational cost. Distillation defines a smearing function, 

-^vccs 

a,(t)= j;F(AW)eF(t)4'=)t(t), (3.4) 
fc=i 

where in this work A*^^-* and ^^^^ are the eigenvalues and eigenvectors of the gauge-covariant 
lattice Laplacian, — V^, evaluated on the background of the spatial gauge-fields of time-slice 
t and sorted by eigenvalue. Here F{\^^^) is some smearing profile and in this work we set 
F = 1. The distillation operator, □, is a projection operator into the subspace spanned by 
these eigenmodes and so = □. Smeared quark fields are constructed by applying this 
distillation operator to each quark field appearing in the interpolating operators, tp = □■0. 
As described in Ref. [22], the outer product structure of □ allows a correlation function 
to be factorised. The computational cost is manageable because inversions are performed 
in the smaller space of distillation vectors and spin, with dimensionality N^ecs^a, where 
Nfj- = 4 is the number of components of a lattice Dirac spinor. Further, this factorisation 
allows the efficient computation of correlation functions with a large basis of operators, 
as well as three-point functions and multi-hadron two-point correlation functions. All the 
quark fields referred to subsequently are smeared using the distillation operator and we 
denote them by ip. 

The simplest interpolating operators are smeared local fermion bilinears of the form 
O = 'ijj{x)T'il>(x), where F is any product of Dirac gamma matrices and for clarity we have 
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Table 2. Gamma matrix naming scheme. 



J A(dim) 

Ml) 

1 Ti(3) 

2 r2(3)es(2) 

3 ri(3)©r2(3) 0^2(1) 

4 ^i(i)©ri(3)©r2(3)©-E(2) 

Table 3. Continuum spins, J < 4, subduced into lattice irreps, A(dim), where dim is the dimension 
of the irrep. 



suppressed spin, flavour and colour indices. However, such operators do not give access to 
J > 2 or even to the complete set of J^*^ quantum numbers for J < 1. In addition, as will 
be discussed later, we require a large basis of operators in each quantum number channel 
in order to reliably and precisely extract the energies of excited states. 

To give us a large basis of operators with a range of spatial structures in each channel, 
we use the derivative-based construction for operators described in Refs. [10, 11]: gauge- 
covariant spatial derivatives are combined with a gamma matrix within a fermion bilinear. 
The operator is of the general form 

i^T^i'^yi:, (3.5) 

where D = is a lattice-discretised gauge-covariant derivative. As discussed in 

Ref. [11], operators with definite J^^ can be constructed using these "forward-backward" 
derivatives. The derivatives and vector-like gamma matrices are first expressed in the 
circular basis so that they transform as J = 1. For creation operators the basis is 



D m=±i = ^—^{D x^iD y), Drn=0 = iDz- (3.6) 



Using the standard SO{3) Clebsch-Gordan coefficients, we can combine any number of 
derivatives with gamma matrices to construct operators with the desired quantum numbers, 
qJM ^ where M is the Jz component. The choice of T and the way in which the derivatives 
have been coupled determine the parity, P, and charge conjugation, C, quantum numbers. 
The notation we use follows that of Ref. [11]: an operator (F x D^j^^Y contains a gamma 
matrix, F, with the naming scheme given in Table 2, and A^" derivatives coupled to spin 
Jd, overall coupled to spin J. 

In lattice QCD calculations the full three-dimensional rotational symmetry of an infi- 
nite volume continuum is reduced to the symmetry group of a cube. Instead of the infinite 
number of irreducible representations labelled by spin, J > 0, we instead have a finite num- 
ber of lattice irreducible representations [irreps); the five single-cover irreps for states at 
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Figure 2. Principal correlator fits according to Eq. (3.9) for a selection of low-lying states in the 
irrcp on the 24^ volume. Plotted are A"(t) • e™"'*"*"^ data and the fits for to = 15 showing 
the central values and one sigma statistical uncertainties; the grey points are not included in the 
fits. From left to right the states have ajm = 0.53726(4), 0.6713(5), 0.676(1) and 0.753(2). In all 
cases, for the to used we obtain reasonable fits with x^/^dof ~ 1- 
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Table 4. The number of operators used in each lattice irrep, A^*^. All combinations of gamma 
matrices and up to three derivatives are included. 



rest are ^i, Ti, T2, E and Note that parity and any flavour quantum numbers such as 
charge conjugation are still good quantum numbers on the lattice. The distribution of the 
various components, Af, of a spin-J state across the lattice irreps is known as subduction 
and in Table 3 we show the pattern of these subductions for J < 4. Because a state with 
J > 1 is split across many lattice irreps and each lattice irrep contains an infinite tower of 
J, it is not straightforward to identify the J of states extracted in lattice computations. 
Refs. [10, 11] presented a method to address this problem and we discuss this further in 
Section 4. 

Refs. [10, 11] also discussed how to construct operators which transform in a definite 
lattice irrep. A, and row, A, from the continuum operators, O"^'^ , 

0Z = T.<m0'''^ ^ (3.7) 

M 

where are subduction coefficients. In this work we include operators built from all 

possible combinations of gamma matrices with up to three derivatives and then subduced 
into lattice irreps. Table 4 lists the numbers of operators in each irrep. 

We use the variational method [25, 26] to extract spectral information from the two- 
point correlation functions with the particular implementation described in detail in Ref. [11] 
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In brief, for each lattice irrep we form a matrix of correlators, Cij{t), where i and j label 
operators in our basis for that irrep. The best extraction (in the variational sense) of the 
energies and overlaps then follows from solving a generalised eigenvalue problem, 

Cij{t)v] = X''it,to)Cij{to)v] , (3.8) 

where an appropriate reference time-slice to is chosen as described in Refs. [7, 11], the 
energies are determined by fitting the dependence of the eigenvalues. A", called the principal 
correlators, on t — Iq, and the overlaps are related to the eigenvectors, v^. 
We fit the principal correlators using the fit form, 

Xn{t) = (1 - ^n)e-'^"(*-*«) + Ane-<(*-*o) , (3.9) 

where the fit parameters are nin, rn'^ and A^- This form allows for a second exponential 
and empirically we find that the contribution of this second exponential decreases rapidly 
as we increase to- Some example fits to the principal correlators are shown in Fig. 2. These 
are plotted with the dominant time-dependence due to state n divided out and so we would 
see a horizontal line at 1 in the case where a single exponential dominates the fit. 

The combination of the variational method, our operator constructions, the distillation 
technique, and these anisotropic lattice ensembles has been shown to be very effective in 
studies of excited light isovector [10, 11] and isoscalar [12] mesons and baryons [14, 15]. In 
the following section we review our method for identifying the continuum spin of extracted 
states and present its first application to charmonium spectroscopy. 

4 Determining the spin of a state 

In principle, the spin of a single-hadron state can be determined by extracting the spectrum 
at various lattices spacings and then extrapolating to the continuum limit. There, where 
full 5*0(3) rotational symmetry is restored, energy degeneracies between different lattice 
irreps should emerge. For example, according to the pattern of subductions (Section 3), 
a spin-2 state would be expected to appear as degenerate energies within the T2 and E 
irreps. However, there are difficulties with this procedure. First, this requires high precision 
calculations with successively finer lattice spacings and so with increasing computation cost; 
this is not practical with currently available resources. Second, the continuum spectrum 
can exhibit degeneracies and this is particularly true in the charmonium system where 
hyperfine splittings are small. The question then arises as to how to identify, without 
infinite statistics, which degeneracies are due to the lattice discretisation and which are 
physical degeneracies. 

An alternative approach, proposed in Refs. [10, 11], is to consider the overlaps, Eq. (3.3), 
using this additional information to determine the spin of extracted states at a single lat- 
tice spacing. Here we will only review the important points. The subduced operators 
constructed in Section 3 respect the symmetries of the lattice. However, we find that each 
operator, O^^^ carries a "memory" of the continuum spin, J, from which it is subduced. If 
our lattice is reasonably close to restoring continuum symmetry at the hadronic scale, we 
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Figure 3. The normalised correlation matrix, (Cy / \/CuCjj) with t/at = 5 in the irrep on the 
24'^ volume. The operators are ordered such that those subduced from spin 1 appear first followed 
by spin 3 and then spin 4. The correlation matrix is observed to be approximately block diagonal. 

can expect that an operator subduced from spin J will only overlap strongly onto those 
states of continuum spin J [27]. 

This is apparent in both the correlation matrix and the overlaps of individual states. 
Fig. 3 shows the correlation matrix for at time-slice 5 where the operators are ordered 
according to the spin from which they are subduced. It can be seen that the correlation 
matrix is approximately block diagonal. In Fig. 4 we show the overlaps of a selection of 
operators onto states in the irreps A^^. The J = operators are colored black, J = 1 are 
red, J = 2 are green, J = 3 are blue and J = 4 are gold. Every state has a clear preference 
to overlap only onto those operators subduced from a single continuum spin. For clarity 
we only show a subset of the operators but this pattern is observed for our full operator 
basis. 

To be more quantitative, we can compare the overlaps obtained in different irreps. 
The continuum operators have definite spin so that (J',M'|C'-^'*^t|o) = zM(5j, j'^M.M' and 
therefore we have(J',M'|oJ{'5|0) = 5;('f' ZM(5j,j/. The value of is common to different 
lattice irreps: because the subduction coefficients form an orthogonal matrix, the spectral 
decomposition will contain terms proportional to Z^'^^*Z^'^^] these do not depend on the A 
subduced into, up to discretisation effects. This suggests that we can identify the spin of a 
state by comparing the Z-values obtained independently in different irreps. For example, 

[J=2l 

a J = 2 state created by a 0\ operator will have the same Z-value in each of the T2 
and E irreps. 

In Fig. 5 we show the Z-values for states conjectured to have J = 2, 3 and 4 in the 
irreps A~~. It can be seen that there is generally good agreement between the Z- values 
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Figure 4. The overlaps, Z, of a selection of operators onto some of the lower-lying states in each 
of the lattice irreps, A , on the 24'^ volume; states are labelled by their mass in lattice units, atm. 
The Z-values are normalized so that the largest value for that operator across all states is equal to 
1. The lighter area at the top of each bar represents the one sigma statistical uncertainty on either 
side of the mean. Different color shading represents different operators. 



extracted in different irreps although there are some deviations from exact equality. The 
various discretisation effects have been discussed previously [11]. 

For the J > 2 states which have components distributed across different irreps, the 
question then arises as to what we should quote as our estimate of the mass of the state. 
The mass obtained from fits to principal correlators in each irrep can be different due 
to discretisation effects and fitting fluctuations. To minimize the fluctuations caused by 
changes in the fitting of the principal correlators we perform a joint fit to the principal 
correlators [11]. There is a common mass, nin, in the fit but we allow the second exponential 
in the principal correlator [Eq. (3.9)] in each irrep to be different and so the fit parameters 
are rrin, {nT-',^} and }• Our final results are determined from such joint fits and we find 
they are generally reasonable, with x^/^dof ~ 1- In Fig- 6, as an example, we present the 
joint fit of the lightest 4~~ state in the four irreps A^~ , T^ , T2 and E . 

In summary, we have demonstrated that the Z-values of our carefully-constructed 
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Figure 5. A selection of Z-values for states conjectured to have J = 2, 3,4 in irreps A on the 
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Figure 6. A simultaneous fit to the four principal correlators of the lightest 4 state on the 24^ 
volume using a common mass, rrin. Plotted are A"(t) • e™"'^*~*°) data and the fit; the grey points 
are not included in the fit. 



subduced operators can be used to identify the continuum spin of the extracted states. 
This was shown in the light meson sector [11] and here we have tested the method in an 
explicit calculation of the charmonium spectrum. When we present our final spin-assigned 
spectra, for J > 2 states we use the results of joint fits to the principal correlators. 

5 Stability of the extracted spectrum 

In this section we present the results of some systematic tests of the stability of the extracted 
spectrum. We consider variations of the reference time-slice, used in the variational 
method, the number of distillation vectors and changes to the spatial clover coefficient, Cc,. 



- 12 - 



' I 



Jt |u 



% a 



§ § § § § 
♦ ❖ « ❖ « 



—1 I i_ 



1 St excited state 



0.53740 

0.53735 

0.53730 

0.53725 

0.53720 

0.53715 

0.53710 
0.678 
0.677 
0.676 
0.675 
; 0.674 
^ 0.673 
0.672 
0.671 
0.670 



Ground state 



2nd and 3rd excited states 



Figure 7. Extracted mass spectrum as a function of tg for the lower-lying states in the irrep 
on the 24:^ volume; the top left pane shows the lightest ten states and the other panes show the 
lightest four states in more detail. Horizontal bands show the masses extracted at to = 15 with the 
width ±1(7 from the mean. For large enough but not too large to the spectrum is seen to be stable 
under changes in Iq. 



5.1 Variational analysis and Iq 

In our implementation of the variational method the choice of a suitable to is guided by 
how well the correlator is reconstructed, as described in Ref. [11]. In this study we find 
that typically to/ at ^ 12 — 15 is appropriate, depending on the irrep being considered. 

Fig. 7 shows the variation with to of the lightest ten extracted masses in the 
irrep on the 24^ ensemble. It can be seen that for large enough but not too large to, the 
extracted spectrum is stable under variations in to- We suggest that this stability is in a 
large part due to the inclusion of a second exponential term in Eq. (3.9). This term can 
'mop up' other states which leak into the principal correlator because our operator basis is 
finite. We find that typically the mass in this second exponential, m' , is larger than that 
of all of the extracted states. As to increases, the contribution of this second exponential 
falls rapidly due to it having a smaller amplitude. A, and a larger m' . 

In Figs. 8 and 9 we show the extracted overlaps as a function of to for a selection of 
the states in the irrep. We observe that the overlaps show somewhat more sensitivity 
to to than the masses. This might be expected because the Z-values are related to the 
eigenvectors of the generalised eigenvalue problem; these are likely to be strongly effected 
by a 'wrong' orthogonality when the number of contributing states is larger than the 
number of operators in the basis. 
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Figure 8. Extracted overlaps, Z, with operator (ai x D^j^^^)jr^^ as a function of for the Hghtest 
three J — 1 states in the irrcp on the 24"^ volume. The Z-values for the 1st and 2nd excited 
states have been arbitrarily scaled by factors of 3 and 1.25 respectively to fit on the plot. Coloured 
bands show fits to an exponential plus a constant or a constant over an appropriate range of to- 
The Z-values are seen to plateau for sufficiently large but can show significant curvature at small 
to. 
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Figure 9. Extracted overlaps, Z, with operator (oi x Dj^^^r^ j=3)t\'^ as a function of to for the 
lightest J = 4 state in the irrep on the 24^^ volume. The horizontal band shows a fit to a 

constant over an appropriate range of to. The Z- value is seen to plateau for large to but shows 
significant curvature at small to. 



In summary, we conclude that our implementation of the variational method gives a 
reliable extraction of the spectrum if to is chosen appropriately. The masses extracted show 
very little variation with to whereas the Z-values can show a more significant dependence 
on to- This was also observed in a study of the spectrum of light isovector mesons [11]. 
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Figure 10. Extracted spectrum for the T-^ irrep on the 16^ volume as a function of the number 
of distiUation vectors, A^vecs, using lower statistics than the main 16^ results. Horizontal bands 
correspond to the masses extracted with A^vecs = 64 and give the 1 sigma range on either side of 
the mean. The spectrum is observed to be stable for A'vccs ^ 48. 



5.2 Varying the number of distillation vectors 

As described in Section 3, in the distillation method there is a choice of how many eigen- 
vectors of the Laplacian, N^ecs, to include in Eq. (3.4). Using fewer distillation vectors 
reduces the computational cost but using too few results in over smearing and limits our 
ability to extract high-spin and excited states [11]; the unsmeared limit corresponds to 
-^vecs = 3(L/as)^ where 3 corresponds to the number of colours. To get the same smearing 
operator on a larger volume requires A^vecs to be scaled by a factor equal to the ratio of 
the spatial volumes [22]. Therefore, to optimise the efficiency of these calculations, it is 
important to use the smallest number of distillation vectors for which the states of interest 
can be extracted reliably. 

In Fig. 10 we show the low-lying states in the irrep on the 16^ volume as a function 
of the number of distillation vectors used. These correlators were calculated using lower 
statistics, 90 gauge-field configurations and 6 time-sources per configuration, compared to 
the main 16'^ results and so the extracted spectrum at Avecs = 64 is not identical to that in 
the full results presented later. The spin of some states for Avecs ^ 32 could not be reliably 
identified and these are not included in the figure. The spectrum is reasonably stable for 
Avecs ^ 48, although the statistical precision reduces slightly with Avecs = 48 compared to 
Avecs = 64. The quality of the spectrum degrades for fewer vectors and this is particularly 
apparent for the excited and high-spin states. A possible explanation for this behaviour 
was given in Ref. [11]. 

In summary, the number of distillation vectors must be sufficiently large in order to 
reliably extract high-spin and excited states, and once the number of vectors is large enough 
the spectrum is stable under changes in the number. These results are consistent with those 
of Ref. [11]. The observation that 48 vectors are sufficient on the 16^ volume suggests that 
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Figure 11. The mass shifts, measured on the 16'^ ensemble, for the fowest-lying states in a selected 
set of charmonium lattice irreps as , the size of the chromomagnetic clover term, is varied from 
tadpole-improved tree- level value of Cg — 1.35 (green) to = 2 (cyan). Masses presented here are 
measured relative to M^^. Experimental data is shown as solid lines; note that experimentally the 
XcO has a significant hadronic width. 



48 (11) = 162 vectors is reasonable for the 24^ volume. 
5.3 0{as) effects and the hyperfine splitting 

An accurate determination of the hyperfine splitting between the ground-state pseudoscalar 
and vector mesons made from two heavy quarks has been a long-standing problem for lat- 
tice investigations. The determination of this energy is known to be delicately sensitive to 
artefacts arising from the discretisation of the Dirac operator on the lattice. The result for 
the J — rjc splitting determined on our 16^ ensemble using tree-level 0{a) improvement 
is 80(2) MeV, consistent with that obtained on the 24^ ensemble and substantially un- 
dershooting the physical value of 117 MeV. This determination neglects the effects of the 
disconnected Wick graphs contributing to the correlation function of charmonium mesons. 
These diagrams are not expected to change the hyperfine splitting substantially [28] due to 
OZI suppression and their absence is unlikely to explain in full the discrepancy noted above. 
We have made some first tests of the effect of disconnected diagrams in our calculation and 
come to the same preliminary conclusion. One other source of systematic uncertainty is 
the use of up and down quarks heavier than their physical values in our lattice action. 
Other lattice calculations [8] suggest this will also result in only a small effect which will 
not explain the discrepancy we see. 

The lattice action for charm quark fields used in this study is C'(a)-improved at 
tree level in tadpole-improved perturbation theory. It would be expected that a non- 
perturbative determination [29] of action parameters would yield a larger value of c^, the 
coefficient of the spatial Pauli-like term, ipaijFijip, appearing in the lattice action, than the 
value Cs = 1.35 we use. Studies show that increasing the value of Cg yields a larger hyperfine 
splitting at finite lattice spacing. To investigate the effect of a larger value of c^, a study of 
some charmonium states using = 2 was carried out on the 16^ volume. To be clear, the 
choice of Cs = 2 is an ad hoc one, chosen from the consideration that a non-perturbative 
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determination is expected to be larger than the tadpole-improved tree-level value which in 
our simulations is c<j = 1.35. 

The lowest lying ^, 1 , (0,1,2)++ and 1 ^ states were computed, corresponding 
to the lightest S and P-wave states in a quark model and the lightest exotic, and Fig. 11 
shows these masses with the ?7c-mass subtracted. Choosing = 2 gives a mass difference 
between the J/tp and r/c of 114(2) MeV, very close to the experimental value and also leads 
to significantly better agreement in the P-wave system. The spin-two Xc2 state is split on 
the lattice across two irreps, T^"*" and E~^~^. For both values of Cg, these two states are 
degenerate within the statistical precision of this study. This is still consistent with the 
explanation that the hyperfine splitting is underestimated due to 0{a) lattice artefacts. 
Since all dimension-five operators consistent with lattice symmetries do not break the 
continuum rotation group, any differences between these lattice states should only arise at 
O(a^) in a Symanzik expansion and so would be expected to be small. 

The lightest charmonium exotic, the 1 state, shows mild dependence on Cg- This 
suggests the energies of such states in our study should be close to their continuum values. 
Note that the results presented later in Section 6 use the tree-level tadpole-improved value 
of Cs, not the boosted version of this test. This computation does however enable us to 
assign an approximate scale of 40 MeV to the size of the systematic uncertainty arising from 
the leading-order 0{as) correction and observe that it does not increase for the higher-lying 
states in the spectrum where this has been tested. However, a more detailed investigation 
of these discretisation effects has not yet been performed and would be required for more 
complete control of this systematic uncertainty. 

6 Results 

In this section we present the results from our lattice calculations. We first give the results 
by lattice irrep and compare the spectra obtained on the two lattice volumes. We then show 
the final results with states labelled by their continuum quantum numbers, J^*^. Following 
on from this, in Section 7 we interpret the results and discuss the hybrid phenomenology 
suggested by them. 

6.1 Results by lattice irrep and volume comparison 

In Fig. 12 we show the results of the variational analysis for each lattice irrep, A^*-^, with the 
spectra obtained from the two volumes side by side. The one sigma statistical uncertainty 
on either side of the mean, as determined from the variational analysis, is indicated by 
the vertical size of the boxes. The colour coding gives the continuum spin determined by 
considering the operator overlaps as described in Section 4. States identified as J = are 
coloured black, J = 1 are red, J = 2 are green, J = 3 are blue and J = 4 are gold. 

By comparing the overlaps between different irreps, we are able to match levels which 
correspond to the different components of a state with definite continuum spin subduced 
into lattice irreps; generally we see no significant discrepancy between the masses deter- 
mined in different irreps. The dense spectrum of states at and above atm ~ 0.67 would 
appear to be impossible to disentangle using solely the extracted masses. An advantage of 
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Figure 12. Extracted spectra by lattice irrep, A-^*^, on the 16"^ (lighter shading) and 24^ (darker 
shading) volumes. The vertical size of each box gives the one sigma statistical uncertainty on either 
side of the mean and the colour coding indicates the continuum spin. 



using anisotropic lattices is that all masses we consider have ajm < 1, which is not easily 
achieved for charmonium on isotropic lattices. The spatial lattice spacing is ~ 0.12 fm 
and this is fine enough for an effective restoration of rotational symmetry at the hadronic 
scale, shown by the fact that we can successfully determine the continuum spin of extracted 
states. 

In Fig. 13 we show the spin-identified spectrum labeled by continuum J^^ with the 
results from the two volumes presented side by side; only the well-determined low-lying 
states are included. The dashed lines indicate the lowest thresholds for open charm decay, 
the non-interacting DD and DgDg levels with the D and Dg masses measured on the 16^ 
ensemble. We observe that there is no significant difference between the two volumes at the 
level of the statistical uncertainties, even above the open-charm thresholds, evidence that 
we are not seeing multi-hadron states in our extracted spectra. We used higher statistics for 
the 24^ calculation compared to the 16^ calculation (see Section 2.1). On the 24^ volume 
we can reliably extract and identify the spin of more states and determine their masses 
with higher statistical precision, and we will therefore quote these as our final results. 

In summary, we are able to interpret our spectra in terms of single mesons of definite 
continuum spin subduced into various lattice irreps. This, along with the lack of any 
significant volume dependence, is evidence that we are not observing multi-hadron states 
in our extracted spectra; the distribution of such states across irreps would have a different 
pattern. Further evidence is provided by observing that we do not extract energy levels 
where we would expect non-interacting two-meson levels to appear. These points have 
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Figure 13. The extracted spin-identified charmonium spectrum labelled by J^*^; the 16'^ (open 
boxes) and 24'^ (filled boxes) volumes agree well. The vertical size of the boxes represents the 
one sigma statistical uncertainty on either side of the mean. The dashed lines indicate the lowest 
non-interacting DD and DgDs levels using the D and Dg masses measured on the 16"^ ensemble. 



been discussed in detail in Ref. [11]. 
6.2 Final spin-identified spectrum 

Our final results, the well-determined states on the 24^ volume, are shown in Fig. 14 along 
with experimental masses taken from PDG summary tables [21]; these results are also 
tabulated in Appendix A. The X(3872) is not shown because its J^'-' (l"*"*" or 2 '') has 
not been determined experimentally. In the plot, the calculated (experimental) mass of 
the r]c has been subtracted from the calculated (experimental) masses in order to reduce 
the systematic error from the tuning of the bare charm quark mass (see Section 2). The 
dashed lines indicate the lowest non- interacting DD and DgDg levels using the D and 
Ds masses measured on the 16'^ ensemble (fine green dashing) and using the experimental 
masses (coarse grey dashing). We note that higher up in the spectrum the states become 
less well determined. This is in part because, although we have a relatively large number of 
operators in each channel, the basis is still restricted in size. In order to better determine 
these states we would need to include operators with different radial structures and, in 
order to determine states with J > 5, different angular structures, for example higher 
numbers of derivatives. 

In the following section we interpret our results and compare with the experimental 
situation, but before that we discuss some other lattice calculations. Ref. [7] presented a 
pioneering quenched lattice QCD calculation of the excited charmonium spectrum which 
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Figure 14. Summary of the charmonium spectrum up to around 4.5 GeV labelled by J . The 
red and green boxes are the masses calculated on the 24^ volume; black lines are experimental 
values from the PDG [21]. We show the calculated (experimental) masses with the calculated 
(experimental) rjc mass subtracted. The vertical size of the boxes represents the one sigma statistical 
uncertainty on either side of the mean. The dashed lines indicate the lowest non-interacting DD 
and DgDs levels using the D and Ds masses measured on the 16'^ ensemble (fine green dashing) 
and using the experimental masses (coarse grey dashing). 

was interpreted further in Ref. [30]. In comparison, apart from being a dynamical calcu- 
lation, the current work uses a larger basis of operators and these operators have been 
constructed so as to facilitate the identification of the continuum spin of the extracted 
states. This means that we can precisely extract a much larger number of excited states 
and states with high spin, and reliably identify their continuum J^^ . 

Results from an Nf = 2 calculation (dynamical up and down quarks but quenched 
strange quarks) of excited charmonia are presented in Ref. [9] using isotropic lattices with 
a few different lattice spacings and m-,^. They also consider mixing with light mesons in 
the pseudoscalar channel and mixing with some low-lying multi-mesons states in a few 
channels. In addition, preliminary results from dynamical (A'^^ = 2 + 1) calculations of 
excited charmonium spectra have been presented in Ref. [8]. A range of pion masses and 
lattice spacings are considered and chiral and continuum extrapolations are attempted. 
However, both these references use a smaller operator basis compared to ours making a 
robust identification of the continuum spin of the extracted states difficult. Therefore a 
limited number of states can be reliably extracted and in particular the lightest exotic 
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Table 5. Supermultiplets for quark- antiquark pairs with spin S and relative orbital angular mo- 

p c 

mentum L. Also shown are some hybrid supermultiplets, discussed in Section 7.2, where Jg" " are 
the quantum numbers of the gluonic excitation; exotic J^'^ are shown in bold. 

hybrid (1 '') is difficult to unambiguously disentangled from a non-exotic spin four state. 
7 Interpretation of the results 

In this section we give some interpretation of our results, compare with the experimental 
situation and then discuss hybrid meson phenomenology in more detail. Many of the states 
with non-exotic J^'~' in Fig. 14 appear to follow the n?^~^^Lj pattern predicted by quark 
potential models, where 5 is the spin of the quark-antiquark pair, L is the relative orbital 
angular momentum, J is the total spin of the meson and n is the radial quantum number. 
The quantum numbers of the members of various quark-antiquark supermultiplets are listed 
in Table 5. The '^^'^^Lj assignments for our extracted states are determined by considering 
the operator-state overlaps [30] . 

In the left panel of Fig. 14 we show the negative parity states. We observe the ground- 
state S-wave pair [0 '", 1 ] and at M — M^^ ~ 700 MeV the first excitation; a second 
excitation is observed at M — M^^ ~ 1150 MeV. There is a complete D-wave set [(1, 2, 3) , 
2 ^] just above the DD threshold. In the region above 1200 MeV there is a complete set 
of D-wave excitations and parts of a G-wave indicated by the presence of spin-4 states. In 
addition, there are three states at M - M^^ ~ 1200 to 1300 MeV with J^'^ = (0, 2)-+, 1 — 
which do not appear to fit with the pattern expected by quark models. We note that these 
three states all have a relatively large overlap onto operators proportional to the gluonic 
field-strength tensor, something not observed for the states that fit into quark-antiquark 
supermultiplets. For example, this 1 , the sixth state in the irrep in Fig. 4, is the 
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Figure 15. Overlaps, Zf = (n|Oj|0), for the lightest P-wave (left panel) and D-wave (right panel) 
supermultiplets with operators {tt, /o}nr x -DjLi and {tt, p}nr, x D^fL2 respectively. 



only state shown which has a large overlap onto the (vr x DjJ^j^) operator, an operator 
proportional to the chromomagnetic part of the field-strength tensor. Following Ref. [11] 
we suggest that these 'excess' states can be identified as non-exotic hybrid mesons and we 
return to this in Section 7.2. 

In the middle panel we show the positive parity states. Below DD threshold there is a 
clear P-wave set [(0, 1, 2)++, 1+"]. In the region of M-M^^ ~ 1000 MeV there is a P-wave 
excitation and, slightly higher, an F-wave [(2,3,4)+"'", 3"' ]. The band of states around 
M — Mff^ ~ 1400 MeV probably contains part of the second excitation of the P-wave set 
and several non-exotic hybrids which lie above the lightest negative-parity hybrids; we will 
discuss these hybrid candidates in Section. 7.2. 

The states with exotic J^*^ are presented in the right panel. These cannot consist 
solely of a quark-antiquark pair - more degrees of freedom are necessary. In general, there 
are several possible interpretations such as hybrid mesons where the gluonic field is excited, 
molecular mesons or tetraquarks. The exotic states in our spectrum have significant overlap 
onto operators with non-trivial gluonic structure and this suggests that we can identify 
them as hybrid mesons. Further, the lack of evidence for multi-hadron states in any 
channel, as discussed in Section 6, suggests that these exotics are not multi-hadron states. 
The 1 '"is the lightest with M — M^^ ~ 1200 MeV and nearly degenerate with the three 
non-exotic hybrid candidates [(0, 2) '", 1 ] observed in the negative parity sector. Higher 
in the spectrum there is a pair of states, (0, 2)"' , at M — M^^ ~ 1400 MeV and a second 
2'^ state slightly higher again. In Section 7.2 we discuss the hybrid meson phenomenology 
suggested by our spectra. 

In order to further test the assignment of states into the lightest P and D-wave 
supermultiplets, we follow Ref. [31] and compare operator-state overlaps for different 
states within a supermultiplet. Consider the operators (ttnr x Dj^I^)''^^ with J^^ = 
1+- and (pnr X D^]i^y=^'^^^ with J^^ = (0,1,2)++, where pnr = - 7o) and 
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■""NR = 2'y5{^ ~ 7o)-^ These operators have the structure of a quark-antiquark pair in 
a gauge-covariant version of a P-wave with 5 = (vtnr) or 5 = 1 (pnr)- The operators 
(ttnr X Df^^y=^ with = 2-+ and (/>nr x ^'[fla)-^^^'^'^ with J^^ = (1,2,3)^- have 
the structure of a quark-antiquark pair in D-wave with 5* = or S" = 1. 

In Fig. 15 we show the overlaps of these operators with states suggested to be members 
of the hghtest P-wave (left panel) and D-wave (right panel) supermultiplets, and observe 
that these overlaps are significant. More generally, we find that the overlap onto non- 
relativistic operator combinations [30, 31], e.g. ~ 7i(l — 70) and 75(1 — 70), are much larger 
than the complementary combinations, e.g. ~ 7^(1 + 70) and 75(1 + 70), which would be 
zero in the non-relativistic limit. This might be expected since the charmonium system is 
fairly non-relativistic. 

For states within a given supermultiplet, it is expected that the Z-values for each 
of these operators, projected into the relevant lattice irreps, will be similar because they 
correspond to the same underlying spatial wavefunction with differing internal spin and 
angular momentum couplings. This is seen to be the case in Fig. 15, supporting our 
identification of these states as members of the P and D-wave supermultiplets. We note 
that only qualitative conclusions can be drawn from the Z-values for lattice-regularised 
matrix elements; these require renormalisation to be compared with continuum matrix 
elements and this renormalisation can mix operators. 

In summary, our extracted spectrum features states that follow the n?^~^^Lj pattern 
expected by quark models along with others, having both exotic and non-exotic quantum 
numbers, that do not fit into this pattern. We have argued that these extra states can be 
identified as hybrid mesons and discuss their phenomenology below after we compare our 
results with the experimental situation. 

7.1 Comparison with experimental results 

It can be seen from Fig. 14 that the pattern of our extracted masses is in general agreement 
with the rather limited number of experimentally observed states. We note that states 
above threshold can have large hadronic widths and a conservative approach is to only 
consider our mass values accurate up to the hadronic width [11]. Discrepancies could be 
due to discretisation effects arising from the finite lattice spacing, systematic uncertainties 
within the variational method and because we have unphysically heavy up and down quarks. 
We do not expect the unphysical up and down quarks masses to introduce a large error 
in the determination of lower-lying states [8]. However, a more careful treatment of these 
effects is needed in the vicinity of open-charm thresholds [32]. In addition, we have used 
the tree level value for the spatial clover coefficient and, as discussed in Section 5.3, this 
means that we underestimate the S-wave hyperfine splitting. Further sources of systematic 
uncertainty arise from not including disconnected contributions, though these are expected 
to give only a small contribution in charmonium. We have shown in Section 6.1 that we 
do not see any finite volume effects at our level of statistical precision. 

^these are for creation operators and we correct a typo appearing in the definition ttnr in Ref. [3f] 
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In the 1 channel we identify the iIj{2S) with our excited ^Si and the ^/^(3770) with 
our ground state ^-Di. For these extracted states, we find that the admixture of '^Di in the 
dominantly ^Si state is smah, as is the admixture of ^Si in the ^Di state. We could not 
extract a significant signal for the S — D mixing angle. 

Higher up in this channel there are two states in the PDG review, '0(4040) and ■0(4160), 
whose robustness is not beyond question. These states appear in the mass region where 
we have the second excitation of the ^5*1 with a relatively large uncertainty on the mass. 
In light of the observations in Section 5.2, a larger number of distillation vectors may be 
required in order to extract this state more precisely. Alternatively, we may need to add 
operators with more radial structures to our basis, for example by including operators with 
additional spatial derivatives or different smearing profiles. 

The y(4260) at M - M^^ ~ 1300 MeV with J^'^ = f— is one of the states that 
appears to be supernumerary compared to the pattern expected by quark models. A 
possible interpretation is as a non-exotic hybrid meson [33-35]. The mass of the 1 
hybrid candidate in our calculation agrees well with the mass of the y(4260), supporting the 
hybrid interpretation. However, solely comparing masses does not allow strong conclusions 
to be drawn. We find conventional 1 charmonia in the same mass region and from our 
calculation we also can not exclude the possibility that the y(4260) is a multi-meson state 
or a tetraquark. 

We find that the 1 hybrid candidate has a significant overlap onto an operator which 
has the structure of a colour-octet quark-antiquark pair with S = in S-wave coupled to 
the gluonic field, (vr x -DjIi)-^=^. This is in contrast to the conventional 1 mesons which 
have S = 1 and is interesting because the y(4260) decays to vr+vr^J/i/; [21]. The folklore 
was that such decays should conserve the spin of the heavy quarks and so, because the 
J/^p has S = 1, this would imply that the y(4260) also has 5 = 1, at odds with our 5 = 
hybrid. However, the observation by CLEO [36] of a cross section for e~^e~ — )■ 7r+7r~/ic 
(the he has S = 0) comparable to that of e^e'^ —J- tt^tt^ J/ip at a center of mass energy 
of 4170 MeV suggests that this folklore may not be correct. Even so, the y(4260) being 
produced in e~^e~ annihilation suggests that it has some admixture. In order to test 
the interpretation of the y(4260) and draw firmer conclusions about the nature of these 
vector mesons, the unsmeared decay constants ~ (l |V'7jV'| O) could be calculated for all 
the vector states we extract; this is beyond the scope of the present work. In Section 7.2 
we discuss the more general hybrid phenomenology suggested by our results. 

The X(3872) at M-M^^ ^ 890 MeV, very close to DD* threshold, is the 'unexpected' 
state that has been confirmed by the largest number of experiments and perhaps is one of 
the most enigmatic; its structure remains unclear and is the subject of much discussion. 
Possible J^^ values for the X(3872) are 2 and l"*""*". In our results, the D-wave 2 
state has M — M^^ around 30 MeV below the X(3872), while the first excitation of the 
P-wave l^"*" is around 110 MeV above the X(3872). Another possible interpretation for 
the X(3872) is as a molecular state of D and D* mesons. As we have argued, we do 
not appear to see multi-hadron states in our extracted spectrum and so we would not 
expect to observe such a molecular state in this calculation. To pin down the nature of the 
X(3872) we plan to supplement our basis with operators constructed to overlap well with 
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Figure 16. Charmonium spectrum up to around 4.5 GeV showing only J channels in which we 
identify candidates for hybrid mesons. Red (dark blue) boxes are states suggested to be members 
of the lightest (first excited) hybrid supermultiplet as described in the text and green boxes are 
other states, all calculated on the volume. As in Fig. 14, black lines are experimental values 
and the dashed lines indicate the lowest non- interacting DD and DgOg levels. 

multi-meson states [24, 37]. 

7.2 Exotic mesons and hybrid phenomenology 

In Fig. 16 we show the charmonium spectrum for the subset of J^'^ channels in which, 
by considering operator-state overlaps, we identify candidate hybrid mesons. A state is 
suggested to be dominantly hybrid in character if it has a relatively large overlap onto an 
operator proportional to the commutator of two covariant derivatives, the field-strength 
tensor. We note that within QCD non-exotic hybrids can mix with conventional charmonia. 
We find that the lightest exotic meson has J^*^ = 1 and is nearly degenerate with the 
three states observed in the negative parity sector suggested to be non-exotic hybrids, 
(0,2) ^, 1 . Higher in mass there is a pair of states, (0,2)^ , and a second 2"' state 
slightly above this. Not shown on the figures, an excited 1 ^ appears at around 4.6 GeV, 
there is an exotic 3 ^ state at around 4.8 GeV and the lightest exotic does not appear 
until above 5 GeV. 

The observation that there are four hybrid candidates nearly degenerate with J^'-^ = 
(0, 1,2) ^, 1 , coloured red in Fig. 16, is interesting. This is the pattern of states pre- 
dicted to form the lightest hybrid supermultiplet in the bag model [38, 39] and the P-wave 
quasiparticle gluon approach [40], or more generally where a quark-antiquark pair in S- 
wave is coupled to a 1"' chromomagnetic gluonic excitation as shown Table 5. This is not 
the pattern expected in the fiux-tube model [41] or with an S-wave quasigluon. In addition, 
the observation of two 2^ states, with one only slightly heavier than the other, appears 
to rule out the fiux-tube model which does not predict two such states so close in mass. 
The pattern of J^^ of the lightest hybrids is the same as that observed in light meson sec- 
tor [11, 31]. They appear at a mass scale of 1.2 — 1.3 GeV above the lightest conventional 
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charmonia. This suggests that the energy difference between the first gluonic excitation 
and the ground state in charmonium is comparable to that in the hght meson [31] and 
baryon [15] sectors. 

To explore this hypothesis of a lightest hybrid multiplet further, we follow Ref. [31] and 
consider in more detail operator-state overlaps. The operators (/Onr x -Djl^)''^'^'^'^ with 
JPC = (0, 1,2)-+ and (ttnr x Df^^y=^ with J^^ = 1— are discussed in that reference. 
These operators have the structure of colour-octet quark-antiquark pair in S-wave with 
S = 1 (pnr) or 5 = (vtnr), coupled to a non-trivial chromomagnetic gluonic field with 

PC 4_ 

Jg^ " = ^ where J^, Pg and Cg refer to the quantum numbers of gluonic excitation. 
Fig. 17 shows that the four states suggested to form the lightest hybrid supermultiplet have 
considerable overlap onto operators with this structure. 

For states within a given supermultiplet, it is expected that the Z- values for each of 
these operators, projected into the relevant lattice irreps, will be similar as discussed above. 
The relevant overlaps presented in Fig. 17 suggest that the four hybrid candidates have 
a common structure, supporting our identification of them as the members of the lightest 
hybrid supermultiplet consisting of S-wave qq coupled to a 1"' gluonic excitation. We 
note again that only qualitative conclusions can be drawn from the Z-values for lattice- 
regularised matrix elements. 

A heavier hybrid supermultiplet composed of a P-wave colour-octet quark-antiquark 
pair coupled to a gluonic field with Jg^ ^ = 1"* would contain states with J^^ = 0"* , 
(l+-)^ (2+-)2, 3+-, 0++, 1++, 2++ as shown in Table 5. In Ref. [31] this was identified 
as the first excited hybrid supermultiplet in the light meson sector. The exotic 0"^ state 
and two 2"' states are observed in our spectrum with small mass splittings and these 
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have relatively large overlap onto operators with the structure of a P-wave qq coupled to 
a gluonic 1"' . In the non-exotic positive-parity sector, we have evidence from similar 
operator overlaps that in the region around M — M^^ ~ 1400 to 1500 MeV there are 0"*"^, 
1++, 2'^^ and 3"' hybrids as well as candidates for three 1^ hybrids. We therefore 
observe candidates for all expected members of this excited hybrid supermultiplet and 
these are coloured dark blue in Fig. 16. To pin these states down more precisely we would 
need to add operators with more derivatives and spatial structures to our basis; this is 
beyond the scope of the current work. 

In summary, we have identified possible lightest and first excited hybrid supermulti- 

plets and suggest that these can be interpreted in terms of qq pairs in, respectively, S-wave 

PC -±- 

and P-wave, coupled to a chromomagnetic gluonic excitation with Jg^ ^ = 1 . The 
appearing at a much higher mass scale (above ~ 5 GeV) may suggest a different type of glu- 
onic excitation, for example a 1~~ S-wave quasigluon, or have some other origin. We note 
that the inclusion of multi-hadron operators in our basis could modify the interpretation 
of states above open-charm threshold. 

8 Summary 

Using distillation and the variational method with a large basis of carefully constructed 
operators, we have successfully computed an extensive spectrum of charmonium mesons 
with high statistical precision and reliably identified the continuum J^'~" of the extracted 
states. The large number of extracted states up to 4.5 GeV, across all PC combinations, 
goes beyond any other dynamical calculation. The spin identification methodology has 
allowed us to identify the continuum spin of the extracted states and this will be crucial 
in order to compare the higher-lying states against observations once this dense spectrum 
has been more fully mapped out by experiments. 

We find no statistically significant volume-dependence and also that the spectrum 
is stable with respect to changes in the details of the variational analysis and variation 
of the number of distillation vectors. Lattice discretisation effects were investigated by 
changing the spatial clover coefficient; however, a more complete determination of this 
systematic uncertainty would require additional lattice spacings. The results presented 
here are for unphysically-heavy up and down quarks corresponding to a pion mass of 
~ 400 MeV at a single lattice spacing; lattice ensembles with lighter up and down quarks 
and further lattice spacings will be considered in future work. In our calculations we include 
only the connected quark line contributions to charmonium correlators. The disconnected 
contributions are expected to be small in charmonium since they are OZI suppressed. 
Nevertheless, in the distillation framework we can efficiently compute these contributions 
and determine their effect on the spectrum. This has already been done in the isoscalar 
light meson spectrum [12] and by applying the same methodology to charmonium we can 
determine the mixings between cc, ss and {uu + dd)j basis states; these calculations are 
ongoing [16]. 

For the first time, we have computed the dynamical spectrum of charmonia with exotic 
quantum numbers (0^ , 1 ^ , 2^ ) up to 4.5 GeV. A determination of these exotic states 
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is particularly interesting because it points to the presence of explicit gluonic degrees of 
freedom. The spectrum of non-exotic states has many features in common with the ■n?^~^^Lj 
pattern expected by quark models along with some states which do not appear to fit into 
such a classification. We have suggested that these extra states can be interpreted as non- 
exotic hybrids and have identified the lightest hybrid supermultiplet consisting of states 
with = (0, 1, 2)-+, 1~, as well as an excited hybrid supermultiplet. This is consistent 
with the pattern observed in lattice calculations in the light meson sector [31] and can 
be interpreted as a colour-octet quark-antiquark pair coupled to a 1^ chromomagnetic 
gluonic excitation. We find that the lightest gluonic excitation has an energy scale of 
1.2 — 1.3 GeV, comparable to that found in the light meson [31] and baryon [15] sectors, 
suggesting common physics. Our results allow an interpretation of the y(4260) as a non- 
exotic vector hybrid meson, but based solely on comparing masses we can not draw a 
definitive conclusion; we have suggested further work to probe its structure. 

We have presented arguments that, as in the study of light isovector mesons [11], we 
see no clear evidence for multi-hadron states. To study such states we plan to enlarge the 
basis of operators to include those with a larger number of fermion fields. In Refs. [24, 37] 
multi-hadron constructions have been presented and the method of Liischer and its exten- 
sions [42-45] was used to compute phase shifts from the denser spectrum of energy levels 
extracted. The efficacy of these constructions was shown in application to tttt scattering 
in isospin-2. By applying this methodology to the charmonium sector we can determine 
the spectrum of multi-hadron states and then use Liischer's method and its extensions to 
compute phase shifts, and so determine the mass and width of resonances. 

As well as the renaissance in charmonium spectroscopy, in recent years there has also 
been renewed interest in the spectroscopy of open-charm D and Dg mesons with unexpected 
states being observed. The methodology presented in this article can be applied to the 
study of these open-charm states and this work is in progress. 

An important goal of the Hadron Spectrum Collaboration is a full understanding of 
charmonium spectroscopy including properties of resonances. Experiments such as BESIII, 
PANDA and those at the LHC will collect vast numbers of cc states to explore the charm 
spectrum above and below threshold, including states with intrinsic gluonic excitations. 
Our successful application of novel methods including distillation and spin-identification 
to determine the single-particle spectrum of charmonium up to 4.5 GeV is a timely contri- 
bution to this effort. 
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Table 6. Summary of the charmonium spectrum calculated on the 24^ ensemble, as presented in 
Fig. 14, with statistical uncertainties shown. The scale has been set using the fi-baryon mass as 
described in Section 2; we note that ~ 400 MeV in these simulations. The calculated i]c mass 
has been subtracted in order to reduce the systematic error from the tuning of the bare charm 
quark mass. The masses of states with J > 2 are from joint fits to principal correlators as described 
in Section 4. 
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A Table of results 

In Table 6 we tabulate the masses presented in Fig. 14. The calculated mass of the rjc has 
been subtracted in order to reduce the systematic error from the tuning of the bare charm 
quark mass (see Section 2). 
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